Métodos Funcionales en Aprendizaje Automático

Homework 03 - Diffusion Maps

María Barroso y Gloria del Valle

---

Initial Configuration

This cell defines the configuration of Jupyter Notebooks.

%%html
<style>
    .qst {background-color: #b1cee3; padding:10px; border-radius: 5px; border: solid 2px #5D8AA8;}
    .qst:before {font-weight: bold; content:"Exercise"; display: block; margin: 0px 10px 10px 10px;}
    h1, h2, h3 {color: #5D8AA8;}
    .text_cell_render p {text-align: justify; text-justify: inter-word;}
</style>

%matplotlib inline
%load_ext autoreload
%autoreload 2

This cell imports the packages to be used.

import numpy as np

import matplotlib
import matplotlib.pyplot as plt
from matplotlib import colors

from sklearn.datasets import make_blobs, make_swiss_roll, make_s_curve

from mpl_toolkits.mplot3d import Axes3D

Axes3D

matplotlib.rc("figure", figsize=(15, 5))

seed = 123
my_cmap = plt.cm.Spectral

Introduction

This practical assignment is about implementing the manifold learning method Diffusion Maps, following the scikit-learn template for manifold learning algorithms.

We will design two main functions: one for training the algorithm, obtaining the affinity matrix and the embedded coordinates, and a second one for extending these coordinates for new patterns.

Requirements

The objective of this assignment is to complete the class DM sketched below, which should contain at least the following methods.

__init__(self, sigma, n_components, step=1, alpha=1)

• This is the construction method for the class, with the following parameters:
  • sigma: Kernel parameter for the Gaussian kernel.
  • n_components: dimension of the embedding.
  • step: step in the Markov Chain.
  • alpha: density influence. It should take a value in [0,1].
• This method just store the parameters in the fields of the class, to be used when needed.

fit(self, X, y=None)

• This is the training method, the one that performs the Diffusion Maps algorithm.

• This method should store the affinity matrix and also the eigenvectors to computed the transformation.

fit_transform(self, X, y=None)

• This method returns the embedding coordinates for the training data.
• It should also store the affinity matrix and the coordinates needed to compute the transformation (for example over new, unseen patterns).

transform(self, X)

• This method will obtain approximated coordinates for new, unseen data points.
• It uses for this purpose the Nyström Formula.

Some recommendations

• It should allow to fix all the DM possibilities (different steps, density normalization...). To allow a change in the kernel function could be also a nice idea.
• Implement all the auxiliary functions that you may need, for example, for deciding the best parameter values in each case.
• It could be a nice idea to offer a deterministic output, that do not depend on the sign of the eigenvectors.

Implementation

* Complete the `DM` class below, satisfying the described requirements.

from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.metrics.pairwise import rbf_kernel, laplacian_kernel, pairwise_distances
from sklearn.utils.extmath import svd_flip
from scipy.spatial import distance

class DM(TransformerMixin, BaseEstimator):
    """
        Diffusion Maps.
    """
    def __init__(self,
                 n_components: int = 2 ,
                 kernel_type: str = 'rbf',
                 sigma: (str, float) = 'percentil',
                 step: int = 0, 
                 alpha: float = 0.0,
                 percentil: float = 50,
    ) -> None :
        """
        Parameters
        ----------
        n_components:
            Embedded space dimension.
        kernel_type:
            kernel type to compute the affinity matrix. Can be 'rbf' or 'laplacian' kernel.
        sigma:
            Sigma parameter of rbf or laplacian kernel.
        step:
            Number of steps to advance in the Markov chain.
        alpha:
            Normalization parameter to control the density influence.
        percentil:
            Percentile vale to calculate sigma used when self.sigma == 'percentil'. 
        """
        
        self.n_components = n_components
        self.kernel_type = kernel_type
        self.sigma = sigma
        self.step = step
        self.alpha = alpha
        self.percentil = percentil

    def compute_sigma(self):    
        """ Compute sigma value for the affinity matrix distinguishing each possible method. """
        if self.sigma == 'median':
            self.sigma_ = np.percentile(pairwise_distances(self.X), 50)
        elif self.sigma == 'maximum':
            self.sigma_ = np.percentile(pairwise_distances(self.X), 100)
        elif self.sigma == 'auto':
            self.sigma_ = np.sqrt(self.X.shape[1]/2) # like SVM sklearn gamma = 1/n_features
        elif self.sigma == 'percentil':
            self.sigma_ = np.percentile(pairwise_distances(self.X), self.percentil)
        elif isinstance(self.sigma, float) and self.sigma > 0.0:
            self.sigma_ = self.sigma
        else:
            raise ValueError("Please, enter a correct sigma method: 'median', 'maximun', 'auto' or a direct positive sigma value")
        
    def compute_affinity_matrix(self,
                                X: np.ndarray, 
                                Y: np.ndarray = None
    ) -> np.ndarray:
        """ 
        Compute affinity matrix for X data.
        Parameters
        ----------
        X:
            Data matrix
        y:
            Data matrix [Optional]
        Returns
        -------
        Kernel matrix
        """

        if not hasattr(self, "sigma_"):
            self.compute_sigma()

        if self.kernel_type == 'rbf':
            return rbf_kernel(X, Y, gamma = 1/(2*self.sigma_**2))
        if self.kernel_type == 'laplacian':
            return laplacian_kernel(X, Y, gamma = 1/(2*self.sigma_**2))
        else:
            raise ValueError("Please, enter a correct kernel method: 'rbf' or 'laplacian'.")
            
        
    def fit(self,
            X: np.ndarray, 
            y: np.ndarray = None
    ):
        """
        Compute the embedding vectors for data X
        Parameters
        ----------
        X : 
            Array-like of shape [n_samples, n_features] 
            training set.
        y : 
            Ignored
        Returns
        -------
        self : 
            Returns an instance of self.
        """ 
        
        # Store X for Nystrom approximation
        self.X = X
                
        # Kernel definition with affinity matrix
        K = self.compute_affinity_matrix(X)
        
        # Compute distance matrix with density 
        D = np.sum(K, axis=1)**self.alpha

        # Compute density normalization
        K_alpha = K / D / D.T

        # Compute distance matrix with K normalized
        D_alpha = np.sum(K_alpha, axis=1).reshape(-1,1)

        # Compute stationary distribution
        pi = (D_alpha / np.sum(D_alpha)).reshape(-1,1)

        # P reduced to orthonormal basis by diagonal matrix
            # We consider its conjugated matrix A
        D_sqrt = np.sqrt(D_alpha)
        A = K_alpha / D_sqrt / D_sqrt.T     
        
        # Obtain the SVD decomposition
        U, s, V = np.linalg.svd(A, full_matrices=True)

        # sign correction
        U, _ = svd_flip(U, V)
        
        # Compute final eigenvector of P and eigenvalues
        eigenvectors = U / np.sqrt(pi)
        eigenvalues = s**2
        
        # Discard autovalue 1
        self.eigenvalues = eigenvalues[1: self.n_components+1]
        self.eigenvectors = eigenvectors[:, 1: self.n_components+1]
        
        return self

    def fit_transform(self, 
                      X: np.ndarray, 
                      y: np.ndarray = None
    ) -> np.ndarray:
        """
        Compute the embedding vectors for data X and transform X.
        Parameters
        ----------
        X :
            Array-like of shape [n_samples, n_features]
            training set.
        y : 
            Ignored
        Returns
        -------
        X_red : 
            Array-like, shape (n_samples, n_components)
        """
        
        # fit training set
        self.fit(X)
        
        # Compute the embedding
        self.embedding = (self.eigenvalues**self.step)*self.eigenvectors
        
        return self.embedding
    
    def transform(self, X: np.ndarray) -> np.ndarray:
        """
        Transform X.
        This function is implemented using the Nyström formula.
        Parameters
        ----------
        X : 
            Array-like, shape (n_samples, n_features).
        Returns
        -------
        X_red : 
            Array-like, shape (n_samples, n_components)
        """
        
        # Compute Kernel for the new X data and the X stored before
        K = self.compute_affinity_matrix(X, self.X)
        
        # Compute Distance Diffusion
        D = (K.sum(axis=1)**self.alpha).reshape(-1,1)
        
        # Transition Probability Matrix definition
            # Probability of arriving from i to j in one step
        P = K / D
        
        # Compute nystrom features
        X_red = (P@self.eigenvectors)/self.eigenvalues
        
        return X_red

Experiments with DM

* Obtain some good embedded coordinates for the three training datasets specified below. * Extend its coordinates for the new points.

Definimos un objeto StandardScaler para estandarizar todos los datos que utilizaremos a lo largo del notebook. Así estarán centrados y escalados a 1 y será posible realizar el embedding utilizando Diffusion Maps.

from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()

Dataset 1: two blobs

N = 1000

X, y = make_blobs(
    n_samples=N, n_features=50, centers=2, cluster_std=3.0, random_state=seed
)
X = scaler.fit_transform(X)

y[y != 1] = -1
plt.scatter(X[:, 0], X[:, 1], c=y, cmap=my_cmap)
plt.axis("equal")
plt.show()

N_new = 100
X_new, y_new = make_blobs(
    n_samples=N_new, n_features=50, centers=2, cluster_std=3.0, random_state=seed + 1
)
X_new = scaler.transform(X_new)

Se representan las nuevas coordenadas del embedding para $T=\alpha=0$ con un núcleo gaussiano variando el parámetro $\sigma$.

sigma_parameters = ['auto', 'median', 'maximum']
dm_alpha = [DM(sigma = sigma) for sigma in sigma_parameters]

fig, axes = plt.subplots(1, 3, figsize=(20, 5))
fig.suptitle(r'BLOBS (T=$\alpha$=0)')
for i, dm in enumerate(dm_alpha):
    X_red = dm.fit_transform(X)
    axes[i].scatter(X_red[:, 0], X_red[:, 1], c=y, cmap=my_cmap)
    axes[i].set_title(r'DM RBF ($p$={0}, $\sigma$={1:.2f})'.format(sigma_parameters[i], dm.sigma_))

No observamos ninguna diferencia notable entre los embeddings así que elegimos una representación que utilice la mediana.

dm = DM(sigma = 'median')
X_red = dm.fit_transform(X)

plt.scatter(X_red[:, 0], X_red[:, 1], c=y, cmap=my_cmap)
plt.title(r'DM BRF ($\sigma$={0:.2f})'.format(dm.sigma_))
plt.show()

X_red_new = dm.transform(X_new)

plt.scatter(X_red_new[:, 0], X_red_new[:, 1], c=y_new, cmap=my_cmap)
plt.title(r'DM BRF New Data ($\sigma$={0:.2f})'.format(dm.sigma_))
plt.show()

Observamos en la primera figura el embedding obtenido para los datos de entrenamiento. Vemos que las clases son linealmente separables y encontramos un margen de separación bastante grande.

Para los nuevos datos las clases no son linealmente separables. Si se utiliza un margen muy estrecho, vemos que uno de los puntos azules pertenecería a la clase roja. Sin embargo, podemos concluir que se trata de un buen embedding para datos no observados ya que la separación lineal es casi perfecta.

Dataset 2: the swiss roll

N = 1500
X, color = make_swiss_roll(n_samples=N, random_state=seed)
X = scaler.fit_transform(X)

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=color, cmap=my_cmap)
plt.show()

N_new = 100
X_new, color_new = make_swiss_roll(n_samples=N_new, random_state=seed + 1)
X_new = scaler.transform(X_new)

En primer lugar, comprobaremos que los embeddings obtenidos por nuestra clase implementada se corresponden con los mismos que devuelve SpectralEmbedding cuando $T=\alpha=0$ para el núcleo gaussiano.

from sklearn.manifold import SpectralEmbedding

dm = DM(sigma=1.0)
X_red = dm.fit_transform(X)

plt.scatter(X_red[:, 0], X_red[:, 1], c=color, cmap=my_cmap)
plt.title(r'DM RBF ($\sigma$={0})'.format(dm.sigma_))
plt.show()

se = SpectralEmbedding(affinity='rbf', gamma=0.5)
X_red_sk = se.fit_transform(X)

plt.scatter(X_red_sk[:, 0], X_red_sk[:, 1], c=color, cmap=my_cmap)
plt.title(r'Spectral Embedding RBF ($\sigma$={0})'.format(dm.sigma_))
plt.show()

Vemos que coinciden, salgo por el signo y la escala.

A continuación, se representa el embedding de Diffusion Maps para $T=\alpha=0$ con un núcleo gaussiano variando el parámetro $\sigma$.

sigma_parameters = ['auto', 'median', 'maximum']
dm_sigma = [DM(sigma = sigma) for sigma in sigma_parameters]

fig, axes = plt.subplots(1, 3, figsize=(20, 5))
fig.suptitle(r'SWISS ROLL (T=$\alpha$=0)')
for i, dm in enumerate(dm_sigma):
    X_red = dm.fit_transform(X)
    axes[i].scatter(X_red[:, 0], X_red[:, 1], c=color, cmap=my_cmap)
    axes[i].set_title(r'DM RBF ($p$={0}, $\sigma$={1:.2f})'.format(sigma_parameters[i], dm.sigma_))

No se obtienen buenos embeddings debido a la elección del parámetro sigma. Esto se debe a que el parámetro $\sigma$ es demasiado grande y estamos teniendo en cuenta la información de vecinos pertenecientes a otras clases. En este conjunto de datos, la estructura geométrica es más complicada y necesitamos que el embedding utilice información muy local.

Por ese motivo, centraremos la elección de $\sigma$ variando el percentil seleccionado para distintos parámetros de densidad $\alpha=\{0.0, 0.3, 0.5\}$.

alpha_parameters = [0.0, 0.3, 0.5]
percentil_parameters = [0.2, 0.3, 0.5, 1.0, 1.5]

for alpha in alpha_parameters:
    dm_percentil = [DM(percentil = percentil, alpha = alpha) for percentil in percentil_parameters]

    fig, axes = plt.subplots(1, 5, figsize=(25, 5))
    fig.suptitle(r'SWISS ROLL ($\alpha$={0})'.format(alpha))
    for i, dm in enumerate(dm_percentil):
        X_red = dm.fit_transform(X)
        axes[i].scatter(X_red[:, 0], X_red[:, 1], c=color, cmap=my_cmap)
        axes[i].set_title(r'DM RBF ($p$={0}, $\sigma$={1:.2f})'.format(percentil_parameters[i], dm.sigma_))

Observando las distintas configuraciones representadas podemos concluir que la mejor representación de los datos, aunque sin conseguir desenrollar el swiss roll, se obtiene al utilizar

• {$\alpha$ = 0.5, $\sigma$ = 0.244 $\sim$ percentil = $0.5$}

El resto de representación o bien agrupan todos los datos cerca de cero o bien no consiguen eliminar la estructura enrollada de los datos.

A continuación representaremos los embeddings obtenidos al utilizar el núcleo laplaciano realizando el mismo experimento.

alpha_parameters = [0.0, 0.3, 0.5]
percentil_parameters = [0.2, 0.3, 0.5, 1.0, 1.5]

for alpha in alpha_parameters:
    dm_percentil = [DM(percentil = percentil, alpha = alpha, kernel_type = 'laplacian') for percentil in percentil_parameters]

    fig, axes = plt.subplots(1, 5, figsize=(25, 5))
    fig.suptitle(r'SWISS ROLL ($\alpha$={0})'.format(alpha))
    for i, dm in enumerate(dm_percentil):
        X_red = dm.fit_transform(X)
        axes[i].scatter(X_red[:, 0], X_red[:, 1], c=color, cmap=my_cmap)
        axes[i].set_title(r'DM Laplacian ($p$={0}, $\sigma$={1:.2f})'.format(percentil_parameters[i], dm.sigma_))

Al utilizar el núcleo Laplaciano las conclusiones que podemos extraer cambian. En este caso, la mejor configuración que obtenemos, aunque no consigue desenrollar los datos, consigue proyectarlos en el eje x al utilizar $\alpha=0.0$ y $\sigma=0.13$, seleccionado a partir del percentil $p=0.2$. Es decir, que al utilizar un núcleo Laplaciano utilizar información más local consigue mejores resultados.

Ahora, sobre las dos mejores configuración de parámetros obtenidas

• {kernel gaussiano, $\alpha=0.5$ $\sim$ percentil = $0.5$}
• {kernel laplaciano, $\alpha=0.0$ $\sim$ percentil = $0.2$}

ilustraremos el embedding y la aproximación de Nystrom para nuevos datos.

#### Best configuration RBF
dm = DM(percentil = 0.5, alpha = 0.5)
X_red = dm.fit_transform(X)

plt.scatter(X_red[:, 0], X_red[:, 1], c=color, cmap=my_cmap)
plt.title(r'DM BRF ($\alpha$={0}, $\sigma$={1:.2f})'.format(dm.alpha, dm.sigma_))
plt.show()

X_red_new = dm.transform(X_new)

plt.scatter(X_red_new[:, 0], X_red_new[:, 1], c=color_new, cmap=my_cmap)
plt.title(r'DM BRF New Data ($\alpha$={0}, $\sigma$={1:.2f})'.format(dm.alpha, dm.sigma_))
plt.show()

#### Best configuration laplacian
dm = DM(percentil = 0.2,  kernel_type = 'laplacian', alpha = 0.0)
X_red = dm.fit_transform(X)

plt.scatter(X_red[:, 0], X_red[:, 1], c=color, cmap=my_cmap)
plt.title(r'DM Laplacian ($\alpha$={0}, $\sigma$={1:.2f})'.format(dm.alpha, dm.sigma_))
plt.show()

X_red_new = dm.transform(X_new)

plt.scatter(X_red_new[:, 0], X_red_new[:, 1], c=color_new, cmap=my_cmap)
plt.title(r'DM Laplacian New Data ($\alpha$={0}, $\sigma$={1:.2f})'.format(dm.alpha, dm.sigma_))
plt.show()

Observamos que aunque la aproximación Laplaciana es mucho más exacta, no se comporta bien para datos no observados. En cambio, utilizando el núcleo gaussiano se obtiene una buena representación para datos no observados.

Dataset 3: the S curve

N = 3000
X, color = make_s_curve(N, random_state=seed)
X = scaler.fit_transform(X)

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=color, cmap=my_cmap)
plt.show()

N_new = 100
X_new, color_new = make_s_curve(N_new, random_state=seed + 1)
X_new = scaler.transform(X_new)

En primer lugar, probaremos si se obtiene la misma representación que la obtenida por la clase SpectralEmbedding implementada en sklearn.

dm = DM(sigma=1.0)
X_red = dm.fit_transform(X)

plt.scatter(X_red[:, 0], X_red[:, 1], c=color, cmap=my_cmap)
plt.title(r'DM RBF ($\sigma$={0})'.format(dm.sigma_))
plt.show()

se = SpectralEmbedding(affinity='rbf', gamma=0.5)
X_red_sk = se.fit_transform(X)

plt.scatter(X_red_sk[:, 0], X_red_sk[:, 1], c=color, cmap=my_cmap)
plt.title(r'Spectral Embedding RBF ($\sigma$={0})'.format(dm.sigma_))
plt.show()

Vemos que ambas clases devuelven el mismo embedding para la curva 'S'.

Después, al igual que hicimos con el conjunto de datos anterior, representamos el embedding para $T=0$ con un núcleo gaussiano variando el parámetro percentil en lugar de variando el parámetro sigma ya que interesa mantener la información local de los datos. Además, se prueba para 3 parámetros de densidad $\alpha=\{0, 0.3,0.5\}$.

alpha_parameters = [0.0, 0.3, 0.5]
percentil_parameters = [0.5, 1.0, 1.5, 2.0, 2.5]

for alpha in alpha_parameters:
    dm_percentil = [DM(percentil = percentil, alpha = alpha) for percentil in percentil_parameters]

    fig, axes = plt.subplots(1, 5, figsize=(25, 5))
    fig.suptitle(r'S CURVE ($\alpha$={0})'.format(alpha))
    for i, dm in enumerate(dm_percentil):
        X_red = dm.fit_transform(X)
        axes[i].scatter(X_red[:, 0], X_red[:, 1], c=color, cmap=my_cmap)
        axes[i].set_title(r'DM RBF ($p$={0}, $\sigma$={1:.2f})'.format(percentil_parameters[i], dm.sigma_))

Cuando $\alpha=0$, obtenemos representaciones en forma de 'U' y otras representaciones en las que se mantiene la estructura de 'S' de los datos. Cuando $\alpha=0.3$, ocurre lo mismo aunque los datos están más dispersos. En cambio, cuando $\alpha=0.5$, conseguimos desenrollar la 'S' al aumentar el tamaño de $\sigma$ a 0.48 y obtenemos una banda aplanada.

A continuación, realizaremos la misma exploración para el núcleo laplaciano.

alpha_parameters = [0.0, 0.3, 0.5]
percentil_parameters = [0.5, 1.0, 1.5, 2.0, 2.5]

for alpha in alpha_parameters:
    dm_percentil = [DM(percentil = percentil, alpha = alpha, kernel_type = 'laplacian') for percentil in percentil_parameters]

    fig, axes = plt.subplots(1, 5, figsize=(25, 5))
    fig.suptitle(r'S CURVE ($\alpha$={0})'.format(alpha))
    for i, dm in enumerate(dm_percentil):
        X_red = dm.fit_transform(X)
        axes[i].scatter(X_red[:, 0], X_red[:, 1], c=color, cmap=my_cmap)
        axes[i].set_title(r'DM Laplacian ($p$={0}, $\sigma$={1:.2f})'.format(percentil_parameters[i], dm.sigma_))

Al utilizar el núcleo laplaciano, no hemos encontrado una elección de parámetros adecuada que consiga desenrollar los datos en una banda como la obtenida para el núcleo gaussiano. Aún así, se obtienen buenas representaciones, sobre todo para $\sigma=0.330$ y $\alpha=0.5$.

Finalmente, representaremos el embedding para las dos mejores configuraciones de parámetros de ambos núcleos tanto para los datos embebidos como para nuevos datos utilizando la aproximación de Nystrom.

• {kernel RBF, $\alpha$=0.5, $\sigma$= 0.48 $\sim$ percentil=$2$}
• {kernel Laplaciano, $\alpha$=0.5, $\sigma$= 0.3330 $\sim$ percentil=$1$}

# RBF kernel
dm = DM(percentil = 2, alpha = 0.5)
X_red = dm.fit_transform(X)

plt.scatter(X_red[:, 0], X_red[:, 1], c=color, cmap=my_cmap)
plt.title(r'DM BRF ($\alpha$={0}, $\sigma$={1:.2f})'.format(dm.alpha, dm.sigma_))
plt.show()

X_red_new = dm.transform(X_new)

plt.scatter(X_red_new[:, 0], X_red_new[:, 1], c=color_new, cmap=my_cmap)
plt.title(r'DM BRF New Data ($\alpha$={0}, $\sigma$={1:.2f})'.format(dm.alpha, dm.sigma_))
plt.show()

# Laplacian kernel
dm = DM(percentil = 1, alpha = 0.5, kernel_type = 'laplacian')
X_red = dm.fit_transform(X)

plt.scatter(X_red[:, 0], X_red[:, 1], c=color, cmap=my_cmap)
plt.title(r'DM Laplacian ($\alpha$={0}, $\sigma$={1:.2f})'.format(dm.alpha, dm.sigma_))
plt.show()

X_red_new = dm.transform(X_new)

plt.scatter(X_red_new[:, 0], X_red_new[:, 1], c=color_new, cmap=my_cmap)
plt.title(r'DM Laplacian New Data ($\alpha$={0}, $\sigma$={1:.2f})'.format(dm.alpha, dm.sigma_))
plt.show()

Ambos consiguen una buena representación de los datos no observados, destacando la banda obtenida por el núcleo gaussiano.

1. Do you consider that the embedding obtained is good for the previous datasets? Is it the expected one? Why? 2. How sensible is the method to its hyper-parameters? Specify the best values that you have found and the technique employed for it. * Check how much influence the sample density (you can vary N for each dataset). * Check what happens if the number of steps `steps` grows. 3. In your opinion, which is the main advantage and disadvantage of this method?

Answers:

Para el primer conjunto de datos el embedding obtenido era correcto y esperado ya que las distancias entre clases en el embedding son mayores. Por otro lado, el embedding para datos no observados no es tan preciso y las clases no están tan separadas. Para el conjunto de datos del swiss roll, la mejor representación de los datos de entrenamiento se obtenía para

• {kernel laplaciano, $\alpha=0.0$ $\sim$ percentil = $0.2$}

el cual obtenia una especie de proyección de los datos desenrollados. Sin embargo, vimos que no daba tan buenos resultados para datos no observados. Por eso concluimos que utilizar

• {kernel RBF, $\alpha=0.5$ $\sim$ percentil = $0.5$}

era más conveniente. Finalmente, para la curva 'S', encontramos una configuración de parámetro óptima

• {kernel RBF, $\alpha$=0.5, $\sigma$= 0.48 $\sim$ percentil=$2$},

que conseguía desenrollar la 'S' y obtenía buenas representaciones para nuevos datos.

Respecto a los parámetros de densidad $\alpha$ y número de pasos en el grafo $T$, a continuación se realiza el análisis para determinar como influyen en las representaciones de los dos últimos conjuntos.

La influencia de parámetro $\alpha$ permite controlar la influencia de la densidad de la muestra. Este control es importante ya que cuando el muestreo de la variedad no es uniforme y no se muestrea un número suficiente de datos, no tiene por qué conseguir un buen embedding de los datos. Cuando $\alpha=1$, Diffusion Maps captura la geometría subyacente sin interferir en la densidad de la muestra mientras que cuando $\alpha=0$, la influencia es muy fuerte a menos que la densidad sea uniforme.

A continuación, se muestra un análisis de la influencia del parámetro de densidad $\alpha$ según variamos el número de datos de los dos conjuntos con geometrías complicadas $N = \{250, 500, 1000\}$ utilizando en ambas el sigma elegido por el percentil 1.

N = [250, 500, 1000]
alpha_parameters = [0.0, 0.25, 0.5, 0.75, 1.0]

# Swiss roll
for n in N:
    X, color = make_swiss_roll(n_samples=n, random_state=seed)
    X = scaler.fit_transform(X)
    dm_alpha = [DM(percentil = 1.0, alpha=alpha) for alpha in alpha_parameters]
    fig, axes = plt.subplots(1, 5, figsize=(20, 5))
    fig.suptitle(f'SWISS ROLL (N={n})')
    for i, dm in enumerate(dm_alpha):
        X_red = dm.fit_transform(X)
        axes[i].scatter(X_red[:, 0], X_red[:, 1], c=color, cmap=my_cmap)
        axes[i].set_title(r'DM RBF ($\alpha$={0})'.format(alpha_parameters[i]))

En la imagen anterior podemos observar que cuando no hay suficientes datos, no se obtiene ningún embedding bueno sea cual sea el valor de densidad. Cuando empiezan a aumentar los datos vemos que al utilizar $\alpha$ mayores se consiguen también buenas representaciones.

Veamos que ocurre para el caso de la curva 'S'.

N = [250, 500, 1000]
alpha_parameters = [0.0, 0.25, 0.5, 0.75, 1.0]

# S curve
for n in N:
    X, color = make_s_curve(n, random_state=seed)
    X = scaler.fit_transform(X)
    dm_alpha = [DM(percentil = 1.0, alpha=alpha) for alpha in alpha_parameters]
    fig, axes = plt.subplots(1, 5, figsize=(20, 5))
    fig.suptitle(f'S CURVE (N={n})')
    for i, dm in enumerate(dm_alpha):
        X_red = dm.fit_transform(X)
        axes[i].scatter(X_red[:, 0], X_red[:, 1], c=color, cmap=my_cmap)
        axes[i].set_title(r'DM RBF ($\alpha$={0})'.format(alpha_parameters[i]))

En esta ocasión, con pocos datos se consiguen buenas representaciones, de modo que no se obtienen resultados tan catastróficos cuando aumentar el parámetro $\alpha$.

Finalmente analizaremos la influencia del parámetro $T$. El parámetro $T$ es el número de pasos en el grafo de difusión, y controla la cantidad de componentes relevantes en el embedding:

$$X_{red} = \lambda^{T}\psi.$$

Es decir, a medida que aumenta $T$, la diferencia entre los autovalores aumenta, de manera que las últimas componentes son cada vez más insignificantes.

A continuación, representamos la influencia de $T=\{1, 100, 300\}$ fijando los parámetros $\sigma$ a aquellos obtenidos por las mejores configuraciones, es decir, $p=0.5, \alpha = 0.5$ para el Swiss roll y $p=2, \alpha=0.5$ para la curva S, utilizando el núcleo gaussiano.

step_parameters = [1, 100, 300]

# Swiss roll
N = 300
X, color = make_swiss_roll(n_samples=N, random_state=seed)
X = scaler.fit_transform(X)

dm_step = [DM(percentil=0.5, step=step, alpha = 0.5) for step in step_parameters]
fig, axes = plt.subplots(1, 3, figsize=(20, 5))
fig.suptitle('SWISS ROLL')
for i, dm in enumerate(dm_step):
    X_red = dm.fit_transform(X)
    axes[i].scatter(X_red[:, 0], X_red[:, 1], c=color, cmap=my_cmap)
    axes[i].set_title(f'DM (T={step_parameters[i]})')

# S curve
N = 300
X, color = make_s_curve(N, random_state=seed)
X = scaler.fit_transform(X)

step_parameters = [1, 100, 300]

dm_step = [DM(percentil = 2, alpha=0.5, step=step) for step in step_parameters]
fig, axes = plt.subplots(1, 3, figsize=(20, 5))
fig.suptitle('S CURVE')
for i, dm in enumerate(dm_step):
    X_red = dm.fit_transform(X)
    axes[i].scatter(X_red[:, 0], X_red[:, 1], c=color, cmap=my_cmap)
    axes[i].set_title(f'DM (T={step_parameters[i]})')

Si observamos las imágenes de arriba, vemos que la estructura geométrica del embedding se mentiene y que solo difieren en la escala, siendo cada vez menor.

Conclusiones

En nuestra opinión, la principal ventaja de este método es que se modela utilizando procesos de Markov, lo que permite otorgar una estructura temporal a los datos. Además, el parámetro de densidad permite controlar la influencia de la densidad de las muestras, lo que hace el método más versatil. También presenta todos las ventajas de los algoritmos basados en variedades, que permiten reducir la dimensión de los datos utilizando transformaciones no lineales al asumir que estos viven en variedades en las que se preservan las distancias. Como principal inconveniente, ambos miembros de la pareja coincidimos que es el alto coste computacional que supone tanto volver a computar la matriz de afinidad para datos no observados utilizando la aproximación de Nystrom como en la búsqueda de hiperparámetros.